# this note gives some mix-newton method
function innerproduct(x,y)
    return sum(x.*y)
end
innerproduct([1,2,3],[1,1,1])

function bb1(x0, x1, g0, g1, min_step)
    # this method gives distance in -g1 direction.
    Δx = x1-x0;
    Δg = g1-g0;
    step = innerproduct(Δx,Δx) / innerproduct(Δx,Δg);
    return -step * g1
end
bb1([0,0],[1,1],[0,1],[1,2.3],0.1)

function bb1_modify(x0, x1, g0, g1, min_step)
    # this method starts from xt=x1-t*Δx and moves in direction -gt.
    Δx = x1-x0;
    Δg = g1-g0;
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔxdotΔx = innerproduct(Δx,Δx)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    step = ΔxdotΔx / ΔgdotΔx
    return -t*Δx - step*gt
end
bb1_modify([0,0],[1,1],[0,1],[1,2.3],0.1)

function bb2(x0, x1, g0, g1, min_step)
    # this method gives distance in -g1 direction.
    Δx = x1-x0;
    Δg = g1-g0;
    step = innerproduct(Δx,Δg) / innerproduct(Δg,Δg);
    return -step * g1
end
bb2([0,0],[1,1],[0,1],[1,2],0.1)

#bb2_modify is a quite amazing method.
function bb2_modify(x0, x1, g0, g1, min_step)
    # this method starts from xt=x1-t*Δx and moves in direction -gt.
    Δx, Δg = x1-x0, g1-g0;
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    step = ΔgdotΔx / ΔgdotΔg
    return -t*Δx -step*gt
end
bb2_modify([0,0],[1,1],[0,1],[1,2],0.1)

# it is found that bb_component is not good.
function bb_component(x0, x1, g0, g1, min_step)
    # this method gives distance for each component
    Δx,Δg = x1-x0,g1-g0;
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    len = length(Δx)
    d = Vector(undef,len)
    for j in 1:len
        α = Δx[j] / Δg[j] * gt[j]
        d[j] = -α
    end
    return -t*Δx + d
end
bb_component([0,0],[1,1],[0,1],[1,2],0.1)

function mix_newton_distance(x0,x1,g0,g1,min_step)
    # this method gives distance almost in -g1 direction.
    Δx, Δg = x1-x0, g1-g0
    a,b = innerproduct(-g1,Δx), innerproduct(-g1,Δg);
    r1 = (a/b/2)*(-g1);
    if a*b > 0
        return r1
    elseif a*b < 0
        #println("use conter-direction")
        c = innerproduct(-g0,Δg);
        r_aux = a/(c-b)*Δg;
        return (r_aux - r1)
    else
        println("use min_step")
        return (min_step/norm(g1))*(-g1)
    end
end
mix_newton_distance([0,0],[1,1],[0,1],[1,2],0.1)

function mix_newton_distance_v2(x0,x1,g0,g1,min_step)
    # this method gives distance almost in -gt direction, where gt is orthogonal to Δx=x1-x0
    Δx,Δg = x1-x0,g1-g0
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    g1dotΔg = innerproduct(g1,Δg)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    abs(t)<1/2 ? t = t : t = 0;#println("t is not used")
    gtdotΔx = g1dotΔx - t * ΔgdotΔx
    gtdotΔg = g1dotΔg - t * ΔgdotΔg
    λ = gtdotΔx / gtdotΔg
    if λ > 0
        return -t*Δx - λ*(g1-t*Δg)
    else
        #println("use conter-direction")
        return -t*Δx + λ*(g1-t*Δg) - 2*gtdotΔx / ΔgdotΔg * Δg
    end
end

function sub_newton_distance(x0,x1,g0,g1,min_step)
    # this method gives distance almost in span(-g1,d0) subspace.
    d0 = x1-x0;
    coeff = [innerproduct(g1,g1),
             innerproduct(d0,d0),
             innerproduct(g1,d0),
             innerproduct(g0,g1),
             innerproduct(g0,d0)];
    μ_aux = sqrt(coeff[1]/coeff[2]) / 2 * (coeff[3] - coeff[5]);
    μ1 = coeff[3]/(-coeff[4] + coeff[1] - μ_aux);
    μ2 = coeff[3]/(-coeff[4] + coeff[1] + μ_aux);
    if μ1 > μ2
        μ = μ1
        ν = sqrt(coeff[1]/coeff[2]) / 2 * μ
    else
        μ = μ2
        ν = -sqrt(coeff[1]/coeff[2]) / 2 * μ
    end
    return (-μ * g1 + ν * d0)
end
sub_newton_distance([0,0],[1,1],[0,1],[1,2],0.1)

function newton_like_distance(x0,x1,g0,g1,min_step)
    #this method gives distance in span(Δx,g0,g1) subspace.
    Δx,Δg = x1-x0, g1-g0
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    g1dotΔg = innerproduct(g1,Δg)
    ΔgdotΔg = innerproduct(Δg,Δg)
    g1dotg1 = innerproduct(g1,g1)
    t = g1dotΔx / ΔgdotΔx #gt orthogonal to Δx
    τ = g1dotΔg / ΔgdotΔg #gτ orthogonal to Δg
    gt = g1 - t*Δg
    gt_norm = sqrt(innerproduct(gt,gt))
    Δg_norm = sqrt(ΔgdotΔg)
    #Δg_norm_p = sqrt((ΔgdotΔx / Δx_norm)^2 + (Δgdotgt / gt_norm)^2)
    #Δx_norm = sqrt(innerproduct(Δx,Δx))
    #Δgdotgt = innerproduct(Δg,gt)
    δ = gt_norm / Δg_norm
    if t==τ
        return -t*Δx - min_step*gt
    else
        ξ = t + sign(t-τ)*δ #find the point near smaller axes
        λ = (ξ-t)*ΔgdotΔx / (g1dotΔg - ξ* ΔgdotΔg)
        return -ξ*Δx + λ*(g1 - ξ*Δg)
    end
end
newton_like_distance([0,0],[1,1],[0,1],[1,2],0.1)

function newton_like_distance_improve(x0,x1,g0,g1,min_step)
    #this method gives distance in span(Δx,g0,g1) subspace.
    Δx,Δg = x1-x0, g1-g0
    g1dotg1 = innerproduct(g1,g1)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    g1dotΔg = innerproduct(g1,Δg)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx #gt orthogonal to Δx
    τ = g1dotΔg / ΔgdotΔg #gτ orthogonal to Δg
    if t == τ
        println("use min_step")
        return -t*Δx - min_step*(g1-t*Δg)
    else
        gt_norm = sqrt(g1dotg1 - 2*t*g1dotΔg + t^2 * ΔgdotΔg)
        gtau_norm = sqrt(g1dotg1 - 2*τ*g1dotΔg + τ^2 * ΔgdotΔg)
        ϕ = acos(gtau_norm / gt_norm)
        γ = (1+2*tan(ϕ)*(tan(ϕ)-sec(ϕ))) * 0.8 # *0.5
        Δ = 1 - 4*γ/(1-γ)^2 * tan(ϕ)^2
        ξq = -1 + 2/(1-γ)/(1+sqrt(Δ)) #ξq = -1 + 2/(1-γ)/(1-sqrt(Δ))(this is a wrong choice)
        δx,δg = (τ-t)*Δx, (τ-t)*Δg
        gq =  g1+ (-t+ ξq*(τ-t))*Δg
        gqdotΔg = g1dotΔg + (-t+ ξq*(τ-t))*ΔgdotΔg
        return (-t + ξq*(τ-t))*Δx -
            ξq*(τ-t)*ΔgdotΔx/gqdotΔg * (g1+ (-t+ ξq*(τ-t))*Δg)
    end
end
newton_like_distance_improve([0,0],[1,1],[0,1],[1,2],0.1)

function newton_like_distance_improve_v2(x0,x1,g0,g1,min_step)
    #this method gives distance in span(Δx,g0,g1) subspace.
    Δx,Δg = x1-x0, g1-g0
    g1dotg1 = innerproduct(g1,g1)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    g1dotΔg = innerproduct(g1,Δg)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx #gt orthogonal to Δx
    τ = g1dotΔg / ΔgdotΔg #gτ orthogonal to Δg
    if t == τ
        println("use min_step")
        return -t*Δx - min_step*(g1-t*Δg)
    else
        gtau = g1-τ*Δg
        δx = (τ-t)*Δx
        δxdotgtau = innerproduct(δx,gtau)
        if δxdotgtau > 0 #saddle point
            println("close to saddle, and escape!")
            return xxxxxxxxxxxxxxxxxxxxxx
        else
            gt_norm = sqrt(g1dotg1 - 2*t*g1dotΔg + t^2 * ΔgdotΔg)
            gtau_norm = sqrt(g1dotg1 - 2*τ*g1dotΔg + τ^2 * ΔgdotΔg)
            ϕ = acos(gtau_norm / gt_norm)
            γ = (1+2*tan(ϕ)*(tan(ϕ)-sec(ϕ))) * 0.8 # *0.5
            ξp = -1/(1-γ)
            gpdotΔg = g1dotΔg + (-t+ ξp*(τ-t))*ΔgdotΔg
            return (-t + ξp*(τ-t))*Δx -
                ξp*(τ-t)*ΔgdotΔx/gpdotΔg * (g1+ (-t+ ξp*(τ-t))*Δg)
        end
    end
end
newton_like_distance_improve_v2([0,0],[1,1],[0,1],[1,2],0.1)


function newton_like_distance_v2(x0,x1,g0,g1,min_step)
    #this method gives distance in span(Δx,g0,g1) subspace.
    Δx,Δg = x1-x0, g1-g0
    ΔxdotΔx = innerproduct(Δx,Δx)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    return -t*Δx - sqrt(ΔxdotΔx / ΔgdotΔg)*(g1 - t*Δg)
end

function newton_like_distance_v3(x0,x1,g0,g1,min_step)
    Δx,Δg = x1-x0, g1-g0
    ΔxdotΔx = innerproduct(Δx,Δx)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    len = length(x0)
    Δg_subnorm = Vector(undef,len)
    for j in 1:len
        Δgdote2 = Δg[j] - ΔgdotΔx/ΔxdotΔx * Δx[j]
        e2dote2 = 1 - Δx[j]^2 / ΔxdotΔx
        Δg_subnorm[j] = sqrt(ΔgdotΔx ^2 / ΔxdotΔx
                             + Δgdote2 ^2 / e2dote2)
    end
    return -t*Δx - sqrt(ΔxdotΔx)*(g1 - t*Δg)./Δg_subnorm
end
newton_like_distance_v3([0,0],[1,1],[0,1],[1,2],0.1)
newton_like_distance_v3([0,0,0],[1,1,1],[0,1,1],[1,2,1],0.1)

function newton_like_distance_v4(x0,x1,g0,g1,min_step)
    # this method start from xt=x1-t*Δx and move in direction -gt.
    # it is an improvement for newton_like_distance_v2
    Δx,Δg = x1-x0, g1-g0
    ΔxdotΔx = innerproduct(Δx,Δx)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    Δgdotgt = innerproduct(Δg,gt)
    gtdotgt = innerproduct(gt,gt)
    Δg_subnorm = sqrt(ΔgdotΔx ^2 / ΔxdotΔx
                      + Δgdotgt ^2 / gtdotgt)
    return -t*Δx - sqrt(ΔxdotΔx)/Δg_subnorm * gt
end

function newton_like_distance_v5(x0,x1,g0,g1,min_step)
    # this method start from xt=x1-t*Δx and move in direction almost -gt.
    # it is an "improvement" for newton_like_distance
    Δx,Δg = x1-x0, g1-g0
    ΔxdotΔx = innerproduct(Δx,Δx)
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    ΔgdotΔg = innerproduct(Δg,Δg)
    g1dotΔg = innerproduct(g1,Δg)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    gt_norm = sqrt(innerproduct(gt,gt))
    Δgdotgt = innerproduct(Δg,gt)
    gtdotgt = innerproduct(gt,gt)
    Δg_subnorm = sqrt(ΔgdotΔx ^2 / ΔxdotΔx
                      + Δgdotgt ^2 / gtdotgt)
    δ = gt_norm / Δg_subnorm
    α = ((t+δ)*ΔgdotΔx - g1dotΔx) / (g1dotΔg - (t+δ)*ΔgdotΔg)
    β = ((t-δ)*ΔgdotΔx - g1dotΔx) / (g1dotΔg - (t-δ)*ΔgdotΔg)
    return -t*Δx + (α+β)/2 * g1 -(α*(t+δ)+β*(t-δ))/2 * Δg
end


#-- general optimize function
function optimize(f,grad,x_init,
                  min_step, grad_min, max_iters,
                  method_of_distance)
    x0 = x_init;
    g0 = grad(x0);
    f_value = f(x0);
    d0 = - min_step*g0;
    x1 = x0 + d0;
    g1 = grad(x1);
    f_value = f(x1);
    counter =1;
    #println([counter, f_value,x1])
    while sqrt(innerproduct(g1,g1)) > grad_min && counter < max_iters
        d0 = method_of_distance(x0,x1,g0,g1,min_step);
        x0 = x1; g0 = g1;
        x1 = x0+(1-1.5^(-counter))d0; g1 = grad(x1);
        f_value = f(x1);
        counter += 1;
        #println([counter, f_value, x1])
        #println(x1)
    end
    println("iterations = ", counter)
    return [x1, f_value]
end

#=
Above methods do not use the information of the values of f, f0,f1,
Below methods will consider this information.
=#
function newton_like_distance_v6(x0,x1,f0,f1,g0,g1,min_step)
    Δx,Δf,Δg = x1-x0, f1-f0, g1-g0
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    t = g1dotΔx / ΔgdotΔx
    gt = g1 - t*Δg
    gt_norm = sqrt(innerproduct(gt,gt))
    xtdotΔg = innerproduct(x1-t*Δx,Δg)
    gtdotΔg = innerproduct(gt,Δg)
    x1dotg1 = innerproduct(x1,g1)
    x0dotg0 = innerproduct(x0,g0)
    α = (2 * Δf + xtdotΔg - x1dotg1 + x0dotg0) / gtdotΔg
    if α > 0
        return -t * Δx - α * gt
    else
        return -t * Δx - min_step/gt_norm * gt
    end
end
function newton_like_distance_v7(x0,x1,f0,f1,g0,g1,min_step)
    Δx,Δf,Δg = x1-x0, f1-f0, g1-g0
    x1dotΔg = innerproduct(x1,Δg)
    g1dotΔg = innerproduct(g1,Δg)
    x1dotg1 = innerproduct(x1,g1)
    x0dotg0 = innerproduct(x0,g0)
    g1_norm = sqrt(innerproduct(g1,g1))
    α = (2 * Δf + x1dotΔg - x1dotg1 + x0dotg0) / g1dotΔg
    if α > 0
        return -α * g1
    else
        return -min_step/g1_norm * g1
    end
end
function newton_like_distance_v8(x0,x1,f0,f1,g0,g1,min_step)
    #this method gives distance in span(Δx,g0,g1) subspace.
    Δx,Δf,Δg = x1-x0, f1-f0, g1-g0
    g1dotΔx = innerproduct(g1,Δx)
    ΔgdotΔx = innerproduct(Δg,Δx)
    g1dotΔg = innerproduct(g1,Δg)
    ΔgdotΔg = innerproduct(Δg,Δg)
    g1dotg1 = innerproduct(g1,g1)
    t = g1dotΔx / ΔgdotΔx #gt orthogonal to Δx
    τ = g1dotΔg / ΔgdotΔg #gτ orthogonal to Δg
    gt = g1 - t*Δg
    gt_norm = sqrt(innerproduct(gt,gt))
    #Δg_norm = sqrt(ΔgdotΔg)
    Δgdotgt = innerproduct(Δg,gt)
    Δx_norm = sqrt(innerproduct(Δx,Δx))
    Δg_norm_p = sqrt((ΔgdotΔx / Δx_norm)^2 + (Δgdotgt / gt_norm)^2)
    δ = gt_norm / Δg_norm_p
    ξ = t + sign(t-τ)*δ #find the point near smaller axes
    λ = (ξ-t)*ΔgdotΔx / (g1dotΔg - ξ* ΔgdotΔg)
    tor = abs((innerproduct(x1,g1) - innerproduct(x0,g0) - 2*Δf) / innerproduct(x1 - t*Δx, Δg))
    confidence = min(tor,1/tor)
    return confidence*(-ξ*Δx + λ*(g1 - ξ*Δg))
end
function optimize_v2(f,grad,x_init,
                  min_step, grad_min, max_iters,
                  method_of_distance)
    x0 = x_init;
    g0 = grad(x0);
    f0 = f(x0);
    d0 = - min_step*g0;
    x1 = x0 + d0;
    g1 = grad(x1);
    f1 = f(x1);
    counter =1;
    #println([counter, f_value,x1])
    while sqrt(innerproduct(g1,g1)) > grad_min && counter < max_iters
        d0 = method_of_distance(x0,x1,f0,f1,g0,g1,min_step);
        x0,f0,g0 = x1,f1,g1;
        x1 = x0+(1-2.0^(-counter))d0;
        g1 = grad(x1);
        f1 = f(x1);
        counter += 1;
        #println([counter, f_value, x1])
        #println(x1)
    end
    println("iterations = ", counter)
    return [x1, f1, g1]
end



## test-1
function test_f(x)
    G = [1,5,10,20]
    return sum(G .* x .* x) / 2
end
test_f([1,1,1,1])
function test_grad(x)
    G = [1,5,10,20]
    return G.*x
end
test_grad([1,1,1,1])
optimize_v2(test_f,test_grad, [1.0,1.0,1.0,1.0],
            0.058, 0.00001, 10000,
            newton_like_distance_v6)
optimize_v2(test_f,test_grad, [1.0,1.0,1.0,1.0],
            0.058, 0.00001, 10000,
            newton_like_distance_v7)
optimize_v2(test_f,test_grad, [1.0,1.0,1.0,1.0],
            0.058, 0.00001, 10000,
            newton_like_distance_v8)


function test_fun1(method)
    optimize(test_f,test_grad, [1.0,1.0,1.0,1.0],
             0.058, 0.00001, 10000,
             method)
end
# best point is (0,0,0,0)
test_fun1(bb1)#counter=29    -----1st
test_fun1(bb1_modify)#counter=62
test_fun1(bb2)#counter=23    -----1st
test_fun1(bb2_modify)#counter=47   ----- 4th
test_fun1(bb_component)#counter= 9
test_fun1(mix_newton_distance)#counter=82
test_fun1(mix_newton_distance_v2)#counter=104
test_fun1(newton_like_distance)#counter=30  ------1st
test_fun1(newton_like_distance_improve)#counter=29
test_fun1(newton_like_distance_improve_v2)#counter=29
test_fun1(newton_like_distance_v2)#counter=63
test_fun1(newton_like_distance_v3)#counter=43
test_fun1(newton_like_distance_v4)#counter=50
test_fun1(newton_like_distance_v5)#counter=33

## test_2
function test_f2(x)
    3*x[1]^2 + 3*x[2]^2 - x[1]^2 * x[2]
end
function grad_f2(x)
    [6*x[1]-2*x[1]*x[2], 6*x[2] - x[1]^2]
end
function test_fun2(method,init)
    optimize(test_f2,grad_f2,
             init, 0.002, 0.001, 10000,
             method)
end
#这个函数在(3√2,3)处是一个鞍点，极值点在(0,0)
test_fun2(bb1,[1.5,1.5]) #counter=7   ----- 1st
test_fun2(bb1_modify,[1.5,1.5])#counter=7
test_fun2(bb2,[1.5,1.5])#counter=6   ----- 1st
test_fun2(bb2_modify,[1.5,1.5])#counter=7   ----- 1st
test_fun2(bb_component,[1.5,1.5])#counter=9
test_fun2(mix_newton_distance,[1.5,1.5])#counter=16
test_fun2(mix_newton_distance_v2,[1.5,1.5])#counter=14
test_fun2(newton_like_distance,[1.5,1.5])#counter=7  ----- 1st
test_fun2(newton_like_distance_improve,[1.5,1.5])#counter=7
test_fun2(newton_like_distance_improve_v2,[1.5,1.5])#counter=7
test_fun2(newton_like_distance_v2,[1.5,1.5])#counter=8  ----- 1st
test_fun2(newton_like_distance_v3,[1.5,1.5])#counter=8  ----- 1st
test_fun2(newton_like_distance_v4,[1.5,1.5])#counter=8  ----- 1st
test_fun2(newton_like_distance_v5,[1.5,1.5])#counter=9  ----- 1st
#这个函数在(3√2,3)处是一个鞍点，极值点在(0,0)
test_fun2(bb1,[-2,4.0]) #counter=9
test_fun2(bb1_modify,[-2,4.0])#counter=13
test_fun2(bb2,[-2,4.0])#counter=9
test_fun2(bb2_modify,[-2,4.0])#counter=8
test_fun2(bb_component,[-2,4.0])#counter=11
test_fun2(mix_newton_distance,[-2,4.0])#counter=23,漂移至鞍点附近再返回
test_fun2(mix_newton_distance_v2,[-2,4.0])#counter=64,漂移至鞍点附近再返回
test_fun2(newton_like_distance,[-2,4.0])#counter=9
test_fun2(newton_like_distance_improve,[-2,4.0])#counter=10
test_fun2(newton_like_distance_improve_v2,[-2,4.0])#counter=10
test_fun2(newton_like_distance_v2,[-2,4.0])#counter=8
test_fun2(newton_like_distance_v3,[-2,4.0])#counter=8
test_fun2(newton_like_distance_v4,[-2,4.0])#counter=8
test_fun2(newton_like_distance_v5,[-2,4.0])#counter=9
optimize_v2(test_f2,grad_f2,[-2,4.0],0.002, 0.001, 10000,newton_like_distance_v6)
optimize_v2(test_f2,grad_f2,[-2,4.0],0.002, 0.001, 10000,newton_like_distance_v7)
optimize_v2(test_f2,grad_f2,[-2,4.0],0.002, 0.001, 10000,newton_like_distance_v8)
optimize_v2(test_f2,grad_f2,[1.5,1.5],0.002, 0.001, 10000,newton_like_distance_v8)

## test_3
function rosenbrock_n(n)
    function rosen(x)
        sum = 0
        for i in 1:n
            r1 = 10*(x[2i]-x[2i-1]^2)
            r2 = 1 - x[2i-1]
            sum += r1^2 + r2^2
        end
        return sum
    end
end
rosenbrock_n(1)([1.0,2.0])
rosenbrock_n(2)([1.0,2.0,1.0,2.0])
function grad_rosenbrock_n(n)
    function grad_rosen(x)
        grad = Vector(undef,2*n)
        for i in 1:n
            r1 = 10*(x[2i]-x[2i-1]^2)
            r2 = 1 - x[2i-1]
            grad[2i-1] = -40 * r1 * x[2i-1] - 2*r2
            grad[2i] =  20 * r1
        end
        return grad
    end
end
grad_rosenbrock_n(1)([1.0,2.0])
grad_rosenbrock_n(2)([1.0,2.0,1.0,2.0])

function test_fun3(method,n)
    init = Vector(undef,2*n)
    for i in 1:n
        init[2i-1] = -1.2 + 1 * (rand() - 0.5)
        init[2i] = 1.0 + 1 * (rand() - 0.5)
    end
    #println(init)
    optimize(rosenbrock_n(n),grad_rosenbrock_n(n),
             init, 0.002, 0.001, 5000,
             method)
end
function test_fun3_v2(method,n)
    init = Vector(undef,2*n)
    for i in 1:n
        init[2i-1] = -1.2 + 3 * (rand() - 0.5)
        init[2i] = 1.0 + 3 * (rand() - 0.5)
    end
    #println(init)
    optimize_v2(rosenbrock_n(n),grad_rosenbrock_n(n),
                init, 0.002, 0.001, 50000,
                method)
end
test_fun3_v2(newton_like_distance_v6,1000)
test_fun3_v2(newton_like_distance_v8,10000) # iterations ~ 400-600
test_fun3(newton_like_distance,10000)#counter~400~     -----1st
# compare newton_like_distance_v8 and newton_like_distance
function test_fun3_v3(n)
    init = Vector(undef,2*n)
    for i in 1:n
        init[2i-1] = -1.2 + 3 * (rand() - 0.5)
        init[2i] = 1.0 + 3 * (rand() - 0.5)
    end
    #println(init)
    optimize_v2(rosenbrock_n(n),grad_rosenbrock_n(n),
                init, 0.002, 0.001, 50000,
                newton_like_distance_v8)
    optimize(rosenbrock_n(n),grad_rosenbrock_n(n),
             init, 0.002, 0.001, 50000,
             newton_like_distance)
end
test_fun3_v3(10000)


#test# best point is (1,1,...,1)
test_fun3(bb1,10)#counter=fail
test_fun3(bb2,10)#counter=fail
test_fun3(bb_component,10)#counter=fail
test_fun3(mix_newton_distance,100)#counter=1894
test_fun3(mix_newton_distance_v2,100)#fail
test_fun3(newton_like_distance,100)#counter~400-700~     -----1st
test_fun3(newton_like_distance_improve,100)#counter~200-400
test_fun3(newton_like_distance_improve_v2,100)#counter~200-400
test_fun3(newton_like_distance_v2,100)#counter=6000
test_fun3(newton_like_distance_v3,100)#counter=15000
test_fun3(newton_like_distance_v4,100)#counter=7000~9000
test_fun3(newton_like_distance_v5,1000)#counter=5000
test_fun3(bb1_modify,100)#counter=15000
test_fun3(bb2_modify,1)#counter=121~200
test_fun3(bb2_modify,10)#counter=200~400
test_fun3(bb2_modify,100)#counter=500~700
test_fun3(bb2_modify,1000)#counter=500~700
test_fun3(bb2_modify,10000)#counter=500~700
# newton_like_distance_improve and bb2_modify are the best, so compare them carefully.
function test_fun3(n)
    init = Vector(undef,2*n)
    for i in 1:n
        init[2i-1] = -1.2 + 3 * (rand() - 0.5)
        init[2i] = 1.0 + 3 * (rand() - 0.5)
    end
    #println("init = ", init)
    @time println("use bb2_modify method: ",
                  optimize(rosenbrock_n(n),grad_rosenbrock_n(n),
                           init, 0.002, 0.0001, 50000,
                           bb2_modify)[2])
    @time println("use newton_like_distance_improve method: ",
                  optimize(rosenbrock_n(n),grad_rosenbrock_n(n),
                           init, 0.002, 0.0001, 50000,
                           newton_like_distance_improve)[2])
end
test_fun3(1)
test_fun3(10)
test_fun3(100)
test_fun3(1000)
test_fun3(10000)

# newton_like_distance_improve is better than bb2_modify, especially for high-dim problems

#= Comments:
(1)newton_like_distance_improve method is the best
(2)newton_like_distance and bb2_modify method are the second best;
(2)mix_newton_distance and mix_newton_distance_v2 follow;
(3)newton_like_distance_v2 seems like newton_like_distance method, but the efficiency is quite slower than newton_like_distance, which means the mean of intersection-points is necessay, newton_like_distance_v2 approximates an ellipse by a cycle is not a quite good idea, which is why bb method is not good for Rosenbrock function.
=#
